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Abstract 

Quantitative error analysis for simulation of wave propagation in three dimensional 
random media assuming narrow angular scattering are presented for the plane wave 
and spherical wave geometry. This includes the errors resulting from finite grid size, 
finite simulation dimensions, and the separation of the two-dimensional screens along 
the propagation direction. Simple error scalings are determined for power-law spectra 
of the random refractive index of the media. The effects of a finite inner scale are also 
considered. The spatial spectra of the intensity errors are calculated and compared to 
the spatial spectra of intensity. The numerical requirements for a simulation of given 
accuracy are determined for realizations of the field. The numerical requirements for 
accurate estimation of higher moments of the field are less stringent. 


1. Introduction 

To analyze the propagation of waves in random media one normally writes equations for 
the propagation of the statistics of the field (such as the moments). In the case of small- 
angle forward-scattering the differential equations for the propagation of the higher field 
moments are difficult to solve. Consequently one must make various approximations, such as 
the Born and Rytov weak-scattering approximations or the Gaussian-field strong-scattering 
approximation. In fact moments higher than the fourth are so difficult that no solutions are 
known outside of the asymptotic weak and strong approximations. Frequently, then, it is 
necessary to solve the equations by numerical methods. 
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In a numerical simulation, by contrast, one generates a realization of the random medium 
(which has the desired statistics) and calculates the resulting wave field. It is actually a 
numerical experiment rather than a numerical analysis. The desired statistic of the wave 
field can then be estimated by repeating the experiment until the estimation errors are 
satisfactory. The simulation lies conceptually between an experiment and an analysis and 
it offers some unique benefits. 

A single set of simulations will provide all the statistics of the field that are of interest, 
whereas for theoretical and numerical analysis, each statistic must be considered separately. 
This is particularly useful for higher moments or more complex field statistics. It is »»asy to 
change the statistics of the medium and repeat a set of simulations, whereas the calculation 
might require a completely different set of approximations. Furthermore a simulation is an 
excellent check for the sort of errors than can easily occur in a complex analysis. 

A simulation cannot, of course, replace an experiment; but it can help identify the im- 
portant factors in the experiment. For example if the medium is turbulent it is very difficult 
to generate or even to identify periods of uniform homogeneous turbulence. Generally real 
turbulent processes are non-stationary and inhomogeneous in a way which is described as in- 
termittent. It is often very difficult to measure the appropriate statistics of the medium and 
usually impossible to thoroughly sample the region of interest. Thus when an experiment 
disagrees with a theoretical calculation it is difficult to tell if the calculation is at fault or 
if the statistical model of the medium is inadequate. Simulations have been used to resolve 
disagreements between theoretical predictions and experiments 1 . Simulations have recently 
been used to investigate the effects of intermittency on laser propagation by Frehlich 2 . These 
issues are important for predictions of the effects of refractive turbulence on coherent lidar 3,4 . 

Another case of interest occurs when it is impossible to observe the process long enough 
to obtain an ensemble average. This is quite often the case in astronomical observations of 
scattering in the interstellar medium. Here one might wish to know if some feature of the 
observation is due to a peculiarity of the medium (i.e. something deterministic) or if it is 
simply a chance sample from a homogeneous stochastic process. Indeed one may find it very 
difficult to describe the feature in question in terms of moments of the random field. In a 
simulation one can compare the wave field with the realization of the random medium in a 
point by point manner. This approach was used by Coles and Filice 5 to demonstrate the 
effect of refraction by large scale structures on the dynamic spectra of intensity scintillation 
in the solar wind. 

Simulations of wave propagation through a thin phase-screen with only one transverse 
dimension have been widely used in the study of ionospheric propagation, where the statistics 
are believed to be highly anisotropic 6-11 . Simulations of radio propagation through the 
solar-wind use the phase-screen approximation with two transverse dimensions 5,12 However 
in nearly isotropic media, such as the atmosphere, a three-dimensional simulation (two 
transverse dimensions and one propagation dimension) is clearly necessary. These are more 
demanding of computational resources (both memory and processor) and much more care 
must be given to the effect of the computational mesh on the accuracy of the result. These 
issues appear to have been first investigated by Filice 12 ; they have also been discussed 
qualitatively by Martin and Flatte 13 ’ 14 . A review of simulation of wave propagation in 
random media is presented by Martin 15 . For extended random media, Spivack and Uscinski 16 
presented the leading order error scaling for the field and second moment as a function of the 
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separation between the screens A z for simulations in a random medium with one transverse 
dimensions (two-dimensional simulation). 

The purpose of this paper is to present quantitative results for the optimal choice of the 
mesh size for simulation of wave propagation through a three-dimensional random medium. 
The propagation problem we consider is small-angle forward-scattering by a phase chang- 
ing medium. The refractive index of the medium is a three-dimensional Gaussian random 
process n(r) which is completely described by its structure function D n { r) or its spatial 
power spectrum $ n (q). The notation of Prokhorov et al 17 will be used throughout. We 
have investigated power-law processes with structure functions of the form D n ( r) = C 2 r“ -1 . 
For such processes the spatial power spectrum is 

$ n (g) = A(a)Clq~ a ~ 2 where A(a) = r(a + l)sin[(a — l)ir/2]/4n 2 . (1) 

Much of the application of this work is to turbulent media with high Reynolds number in 
which an inner-scale is important. To investigate the effect of the inner-scale l Q on simulation 
error we used the universal spectrum for the Kolmogorov exponent a = 5/3 

$ n (q) = 0 . 0330054 C^“ u/ 3 /(< 7 /o) where f(x) = (1 -F aii + a 2 i 2 + a 3 x 3 )exp(— i). (2) 

We used cq =1.4284, a 2 =1.1987, a 3 =0.1414, as estimated from laser scintillation data by 
Frehlich 18 . We have studied three geometries: a plane wave normally incident on a thin slab 
of random medium; a plane wave normally incident on a half-space of uniform homogeneous 
random medium; and a point source embedded in a uniform homogeneous medium. The 
thin slab or screen problem is a canonical case because the extended medium can be treated 
as a stack of slabs. 


2. Forward Scattering Theory 

The choice of computational mesh depends on the important spatial scales in the random 
field in the observation plane. Here we will review the forward scattering theory necessary to 
estimate these scales. If the scattering angle is small it is convenient to write the field with 
respect to a monochromatic plane wave as F( s, z,<) = /(s, z)ex\>[i(kz — u>t)]. Here s is the 
transverse coordinate and k = 2tt/X. The intensity is defined as 7(s, z) = < /( s, z)f"(s, z ) >. 
It is also useful to normalize the field so that the average intensity < /( s, z) > is unity. This 
leads to the parabolic wave equation 

|^(s, z) = V?/(s, z) + ikn{ s, z)f( s, z) , (3) 

where V t is the transverse gradient operator and n(r) is the deviation of the refractive index 
from its mean value. 

Propagation equations for the moments of the field can be derived under the Markov 
approximation 19,20 . This is equivalent to the assumption that the continuous random 
medium can be decomposed into a series of statistically independent thin slabs. These 
slabs must be sufficiently thin that each acting alone would cause only small phase fluctu- 
ations. The conditions for the validity of this approximation have been determined for the 
first moment or average field, the second moment or mutual coherence function, and the 
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fourth moment 20 . The conditions for the validity of the second and fourth moment are very 
weak and for isotropic turbulence essentially equivalent to the small-angle approximation of 
the parabolic equation. We will also make this approximation in the numerical simulation; 
although it is not essential it greatly simplifies the calculations. 

The propagation of the second field moment or mutual coherence function ^(xi,^) = 
< /(xj)/*(xi) > admits an analytical solution for the cases of interest here. If the field is 
wide sense stationary (i.e. is only a function of S12 = Xj— X2) then the Fourier transform of 
T2 is the spatial power spectrum of the field. The bandwidth of this power spectrum sets the 
required sampling interval. The propagation equation for r 2 (xx,X2) is given by Prokhorov 
et al, 17 in terms of D'(s, z), the longitudinal gradient of the wave structure function 

^(X!,X2,2) = ^(V 2 - V 2 )r 2 - ( 4 ) 


D'{ s, z) = 4tt k 2 J J d 2 q[l - cos(q • s)]$ n (q, q 2 = 0, z). (5) 

In the planar geometry the term (V 2 — V 2 )!^ is identically zero and the solution at range 
R can be written 


r 2 (s, R) = exp[— | Z)p(s, i?)] 
where 

D P (s,R)= I* dzD'(s,z). 

Jo 


( 6 ) 

(7) 


A similar solution applies for a spherical wave if Dp{ s, R) is replaced by 

D s (s,R)= [ R dzD'(sz/R,z). 

Jo 


(8) 


The wave structure functions Dp( s) and Ds( s) can also be written as D( s) =< [<^(r) — 
4>(r + s)] 2 > where <f>(r) is the phase on a geometrical path from the source to the receiver. 
So D( s) is a measure of the differential phase variance on a baseline s. The field coherence 
scale s a can be defined by D(s 0 ) = 1. This scale s 0 is inversely related to the rms scattering 
angle because the Fourier transform of T2(s) is the spatial spectrum T^q). A component 
of the spatial spectrum at frequency q corresponds to a plane wave propagating at an angle 
of d — sin _1 (q//:) to the z axis. Thus the width of the angular spectrum can be defined as 
9 0 = l/(ks a ). 

The propagation equation for the fourth moment of the field = < f(xi) /*(x 2) /(x 3) 
/*(£,) > is 


^ = ^(v;-v! + vl-v;) r < -vr ) 

2V = D'{ s 12 ) + D’(s u ) + D'( s 23 ) + D'(s 34 ) - T»'(s 13 ) - D'{ s 24 ). (9) 


The case of greatest interest is when S12 = S34 = 0. Then = < 7(r)/(r + Sj3) > = 
1 + C/(sx 3 ) where (^/(sxa) is the covariance of the intensity fluctuations. This can be solved 
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in weak scattering using the Born approximation and in strong scattering using various 
series expansions (none of which converge quickly). Solutions for pure power-law and for 
Gaussian spectra of refractive index have been summarized by Prokhorov et al. 17 . Solutions 
for power-law spectra with inner scales have been given by Fante 21 for the case of a plane 
wave incident on a homogeneous half space and by Frehlich 22 for the case of a point source 
embedded in a homogeneous medium. The most important characteristics of the intensity 
covariance are outlined below for power-law and related spectra. 

For weak scattering the intensity variance m 2 = C/(0) is less than unity and the shape 
of C/(s) is independent of the level of the refractive index spectrum, but the variance m 2 is 
linearly proportional to this level. In weak scattering the spatial scale of C/(s) is the radius 
of the first Fresnel zone 77 = ( R/k ) 0,5 when the inner scale is less than the Fresnel scale, 
otherwise it is proportional to the inner scale. The Born variance is m 2 = KD(r /) where K 
is a constant of the order of unity which depends weakly on the shape of the spectrum and 
the incident wave (e.g. equations 26 and 38). The Born variance m 2 is a useful measure of 
the strength of scattering even when m 2 >> 1 , that is when m 2 is not a good approximation 
to m 2 . 

In strong scattering a s m 2 — ► 00 , the field approaches a (complex) Gaussian process, the 
mean field approaches zero, and the intensity covariance approaches C/( s) = exp[— Z)(s)], 
i.e. the intensity has unit variance and spatial scale (at 1 /e) of s 0 . The first correction to 
this strong scattering limit adds two terms to C/(s); C T { s) and Cdr( s)* It is conceptually 
helpful to model the intensity as a diffractive process of scale s 0 modulated by a refractive 
process of larger scale s r 23 . If D(s) is isotropic the scale of the refractive process is the 
distance subtended by the scattering angle, that is s r = 0 O R. The product term is Cd r (s) = 
exp[— £)(s)]C' r (s). Thus C T and Cd r carry equal variance, but, as C T is of much larger 
scale, the scale of Cdr is approximately s Q . If D{ s) is anisotropic the refractive process is 
also anisotropic with the same axial ratio, but not, as one might expect, of the form 0 O R. 
Rather, the minor axis of s T is equ? ! to the major axis of 0 O R 24 - 

Three important scales emerge from the analysis of the second and fourth moments. The 
first is the field coherence length s 0 . The second is the scattering disc 0 O R. The third is 
the radius of the first Fresnel zone rj. The three scales are not independent a s rj = s 0 s r . 
Fortunately these scales can be determined from the geometry and the analytical solution 
for the second moment. 


3. The Propagation Calculation 

The propagation of a wave through a random medium according to the parabolic wave equa- 
tion is the limiting case of propagation through a series of discrete phase screens separated 
by free space. A finite number of screens Ns can be used if the phase fluctuations in each 
screen are small. However the simulation of this problem is greatly assisted if these screens 
are statistically independent. This decomposition into independent weak phase screens is 
equivalent to the Markov approximation which is usually made when such problems are 
solved analytically. The effects of correlated screens on simulations with one transverse di- 
mension was shown to be small 16 when the separation between screens was larger than the 
correlation distance of the refractive index fluctuations. Clearly it is possible to simulate a 
medium in which the Markov approximation breaks down and this would be an interesting 
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project. However in the work reported here we have used independent screens. 

The passage of a wave /(s, z) through a thin screen with phase fluctuations <j>(s) at z=0 
can be expressed as /( s, 0 + ) = /(s, 0~) exp[i<^(s)]. The propagation between the screens in 
free-space is expressed as a convolution using the Fresnel integral. For numerical work it 
is efficient to do the convolution in the Fourier transform domain using the Fast Fourier 
Transform (FFT) algorithm (i.e. /*(q, z) = /‘(q, 0) exp[— izq 2 /(2Ar)] where /*(q, z) is the 
Fourier transform of /(s, z)). 

The simulation begins by generating a series of independent phase screens <j> n { s) for 
n = 1, . . . , Ns separated by A z. The field incident on the first screen is deterministic. The 
field incident on a screen /( s, nAz - ) is propagated to the next by: multiplying the incident 
wave by the phase transmittance exp[i<f> n { s)]; Fourier transforming the field; multiplying 
the transform by the free space propagator exp(— iq 2 Az/2A;); and finally retransforming to 
obtain the field incident on the next screen. That is 

/[s, (n + l)Az - ] = Ft~ 1 {Ft[f(s, nAz~) exp(i<?!> r ,(s))] exp[— iq 2 Az/(2fc)]}. (10) 

The field is a bandlimited random process and we can use sampling theory to choose 
the appropriate grid spacing. In general a field of propagating waves is strictly bandlimited 
because q < k. Thus the original field can be reconstructed from samples taken at intervals 
of As < 2ir/(2k) = A/2 without any loss of information. In this case we are considering 
only small angle forward scattering so the field is actually bandlimited much more severely. 
In fact the spatial spectra of interest typically have a form similar to exp(— q 2 /q 2 ). Thus 
there is no significant power for q > q mox , where q max ~ 3q c . Then a sample interval of 
As ~ 7r/q mor ss s 0 will be adequate. 

The grid must be of finite extent and the characteristics of the Discrete Fourier Transform 
(DFT) imposes periodicity of the fields. In choosing this period we are guided by the 
characteristics of the angular spectrum and the intuitive interpretation of the rigorous path- 
integral solution for the fields. In weak scattering (when the Born approximation holds) 
the field in the observation plane is determined primarily by the scattering medium inside 
the first Fresnel zone of rays 25 because these contributions arrive in-phase. The maximum 
extent of the first Fresnel zone is given by r/. In strong scattering the components of the 
scattered field are more randomly phased. In this case, the important region of scattering 
is not the first Fresnel zone but the radius of the scattering disc 9 0 R which is defined by 
the angular spectrum. The scattering disc is larger than the first Fresnel zone if r/ > s Q . 
This implies D(rj) > 1 which is consistent with other definitions of strong scattering based 
on the intensity variance. As it happens the largest scale of the random field (s r ) is the 
radius of the scattering disc so the screen size is also appropriate for the receiving plane. If 
spatial statistics are desired, as is normally the case, then the screen must be large enough 
that each point of interest is at least 0 O R away from the edge of the screen. A similar 
argument based on the spatial intensity spectrum has been used to estimate the required 
screen dimensions 13-15 . We use a square screen of dimensions 4 9 0 R. A rectangular screen 
was also used for comparison. 
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4. Generation of a Random Screen 


The thin slab of random medium in -Az/2 < z < + Az/2 is described by the spatial spectrum 
of refractive index fluctuations $ n (q, z). The two-dimensional phase fluctuations <j>(s,z) are 
given by 


r+Az/2 

<f>(s,z) = k / dz n(s, z). 

J-Az/2 


( 11 ) 


The two-dimensional spatial spectrum of phase fluctuations under the Markov approxima- 
tion is 


**(<!*. q»o) = 2 rfA 2 $„(q,,,q,,q, = 0,z). 


( 12 ) 


The random phase screen is represented as a set of real numbers {4> mn } which are samples 
of a continuous random process <^(s) over a finite range. The process <^(s) is a member of 
an ensemble with power spectral density $,*(q). The sample interval is Ax, A y and the 
simulation dimension is L X) L y . The {<f> mn } are generated in the frequency domain from 
a set of complex random numbers = A pq + iB pq . Here A pq and B pq are independent 
zero-mean Gaussian random variables both with variance 




AxAy 


( 13 ) 


The Gaussian random numbers A pq and B pq are generated from a pair of uniformly dis- 
tributed random numbers using the Box-Muller transformation. Realizations of the phase 
{<£mn} are obtained from the DFT of The real and imaginary parts of <f> mn are un- 

correlated and represent realizations of the desired random process, i.e., the expectation of 
the respective sample spectra are $^(q). The continuous sample spectrum is 


C<t> ^ = TT^Jo Jo y ^exp(-i'q-s)d 2 s | 2 

(AxAy) 2 M N 

4*mn exp[— i(q r mAx + q y nAy)]| 2 . 

q q 

The discrete sample spectrum is then 
C<t>(2np/ L x ,2nq/ L y ) = 


(14) 


(AxA yf |j*j» _ AxAy | fl,| a 


L X Ly 


47r 2 MN 4 t 2 

and substitution of Eq. (13) in (15) demonstrates that < C^(q) > = $^(q) as required 


(15) 


5. Plane Wave Incident on a Thin Phase Screen 

We tested the canonical problem of the plane wave incident on a thin phase screen using 
the universal power-law spectrum given in Eq. (2) for the phase spectrum, i.e., 


<Mq) = 7V“- 2 /(q/ 0 ) 


(16) 
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where T = 2nk 2 A(a)ClAz. -For a pure power-law spectrum (/ 0 = 0) and for 0 < a < 2, 
the strength of scattering mj, the coherence scale s Q and the structure function Djs{ s) are 
given by 

ml = 4nTT(l — a/2) cos(aTr /4)r° /a = K\(a)DTs(^f) (17) 

K\(a) — 2 Q T(1 + q / 2 ) cos ( q 7 t / 4 ) (18) 

Drs(S) = (»/».)“ = • (») 

We have argued that the sample interval (Ax, Ay) should be of the order of the coher- 
ence scale s a and that the screen size should be of the order of the scattering disc 0 O R. 
The coherence scale decreases with increasing strength of scattering and the scattering disc 
increases. Thus the memory and computation time required increase very rapidly with in- 
creasing strength of scattering. It is important to establish the strongest scattering that can 
be simulated with a specified error. 

A reference calculation of was performed using a 2048x2048 point array with Ax = 
Ay = rjf 64 and L y = L x = 64r/. The effect of the finite sample interval or sampling error 
was estimated by resampling the reference screen at Ax = Ay = rjf 32, rjj 16, rj/8 and 
rjj 4, recalculating the field for each resampled screen, and comparing the results of each 
resampled calculation with the reference calculation. For all cases studied the rms error of 
the real part of the field was equal to the rms error of the imaginary part of the field and 
both were half the rms error of the intensity I rma . These tests were done for two exponents 
(a = 1 and 5/3) and three strengths of scattering ( m 2 = 0.1, 1 and 10). 

The error for a pure power-law spectrum should depend only on A s/s c because the 
shape of the spectrum does not depend on the strength of scattering. We tested this by 
plotting I rma versus A s/s 0 in Fig. 1. The error was well described by power-law model, 
[/ rrai = 0.671(As/s o )°' 514 for a = 1 and 7 rm , = 0.382(As/s o )°' 875 for a = 5/3]. Since the 
error is caused by aliasing it should be a steeper function of A s/s Q for steeper spectra. Indeed 
we find (approximately) that 7 rma oc (As/s 0 ) a ^ 2 . We also find that there is essentially no 
change in the rms intensity mdec as As increases, even when the intensity error I rmi is equal 
to mdec- If we write the intensity of the resampled screen as Idee = I + A 7 then 

= Var(I + A/) = m 2 + I? m , + Cm(I , A/) (20) 

where Cov{I , A/) = < (/ — 1)A7 > is the covariance of 7 and A7. Clearly then the 
sampling error A7 is anti-correlated with the intensity 7, i.e. as As increases Idee tends to 
underestimate the intensity peaks and overestimate the intensity nulls. This is exactly the 
behavior one expects of a low-pass filter. 

We also simulated three cases with an exponent of a = 5/3 and an inner scale of l a = 5As. 
These examples were normalized to the pure power-law case with a = 5/3 by setting the 
spectral levels the same at low wavenumber. For this test the spectrum changes with the 
sample interval so I rms should not be expected to scale as A s/s 0 . The results are displayed 
in the bottom panel of Fig. 1. As one would expect, when As > l Q the result is the same as 
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for the pure power-law case. When As < l 0 the error is much less than the pure power-law 
case because the spectrum drops faster at high frequencies and the effects of aliasing are 
greatly reduced. As for the pure power-law spectra, there is essentially no change in the rms 
intensity m* c , even for 7 rm , equal to mj ec . 

The spatial spectrum of intensity fluctuations <I>/(q) for m 2 = 10 and a pure power-law 
phase spectrum (a = 5/3) is shown in Fig. 2. The spectra of the intensity error A I for 
As = rj] 32 and rj/ 16 are overplotted as open symbols. It is interesting to see that the 
sampling errors are almost white, i.e. the errors are independent. These errors are caused 
by aliasing in the angular spectrum, so they are introduced near the Nyquist frequency, 
however after the propagation calculation they have spread almost uniformly throughout 
the intensity spectrum. 

The errors introduced by the finite dimensions of the simulation or windowing errors 
are more difficult to investigate. The random screen generated using the DFT is periodic. 
Originally we created a large reference screen and calculated a reference field. Then we re- 
duced the window by truncating the reference screen. However the errors calculated this way 
increased much more rapidly than we expected (i.e. I rms oc L y /s r ). The problem is clearly 
that the truncation process removes the periodicity in the screen and the discontinuities at 
the edges of the simulation produce addition errors that increase as the observation point 
approaches the edge of the simulation. This additional error would not exist if the screen 
had been generated at the smaller size. 

The windowing error is more easily examined in the transform domain. The propaga- 
tion algorithm conserves power so, with our field normalization, the mean intensity in any 
observing plane is exactly unity. The intensity spectrum is estimated by the DFT so its 
dc component is also exactly unity. However the dc component of the DFT should have 
included the integral over the resolution cell centered on the origin. Thus we have made an 
error in the spectrum which causes an error in the variance, which can be approximated as 

/ " 2 i$ 7 C j bx t'ZSlt j Lj y 

dq r / dqy$/(q r , q^) (21) 

■2Sir/L x J-2Sir/L y 

where 8 is the fraction of the zero-frequency bin that is incorporated into the DC component. 
This error is easily estimated because the integration is confined to very low frequencies 
(q < 2 tt/L) and a low frequency approximation to the intensity spectrum is known 17,26 . 

*/(q) = 4$*(q) sin 2 [q 2 tf/(2*)l exp[-Z>rs(qJ*/*)] ( 22 ) 


If we choose L > n9 0 R, as we have argued intuitively, then the arguments of both the 
exponential and the sin 2 term are small and one can do a power series expansion. The 
integral can then be done term by term. The first term, if L x = L v , is 


Am 2 a(2TT 8r / / L)* a I(a) 
ml x cos(7ra/4)r(l — a/2) 

where 


1(a) = / dx ( dy(x 2 + y 2 ) 1 a ^ 2 
Jo Jo 


(23) 


(24) 
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and /( 5/3) = 0.89627. An example of the windowing error for the pure Kolmogorov spec- 
trum with m 2 = 1.0 is plotted in Fig. 3. Here m 2 is displayed in the upper panel and 
Am = [ml e j — m 2 ] 1 / 2 is displayed in the lower panel, to permit comparison with the sam- 
pling error / rma . The theoretical model [Eqs. (21), (22) and 6 = 0.65] is a dashed line and 
the best-fit power-law function is the solid line. Clearly both fit adequately in the range of 
interest. 

These calculations confirm our intuitive expectation that the sample interval should be 
a fraction of the coherence scale and that the window should be larger than the scattering 
disc. Since the scattering disc is inversely proportional to the coherence scale the number 
of points on the x and y axes must be N oc 6 0 R/s 0 oc 1/s 2 . For a power-law spectrum of 
exponent a the coherence scale s a oc (m 2 ) -1 / a . Therefore the number of points needed is 
N P 2 oc (ml) 4/a . 


6. Plane Wave Incident on a Uniform Medium 

The simulation of an extended medium is performed by collapsing slabs of random media 
with thickness A z into a series of Ns thin weak phase screens as discussed before. The 
sample mesh for each screen must be chosen so Ay ~ 2 rj/y/N and L v ~ 2r/y/N. The 
remaining question is how to choose A z the screen spacing. 

For a pure power-law spectrum ( Iq = 0) and for 0 < a < 2, the strength of scattering 
m 2 , the coherence scale s 0 and the structure function Dp{ s) are given by 22 


ml = K 2 (a)D P (r f ) (25) 

K 2 (a) = 2“ +1 r(l + a/2) cos(a7r/4)/(2 + a) (26) 

Dp( s) = (s/s 0 ) Q = B,(a)k 2 C 2 n Rs a (27) 

B x {a) = 2 1 " £ T(a)r(l - a/2) sin[(a - l)w/2]/r(l + a/2). (28) 


For the Kolmogorov exponent, f?i(5/3) = 2.91438 and K 2 (5/3) = 0.421602. The strength 
of each phase screen is given by Eq. (16) as T = 2irk 2 A(a)ClAz. After the field has 
propagated a distance A z, the perturbation of the field introduced by the phase screen must 
be small. Because the field incident on the phase screen is random, a rigorous estimate of 
the magnitude of the field perturbation is difficult. A simple estimate of the magnitude of 
the field perturbation is the Born variance m 2 [Eq. (17)]. This implies that A z oc 1/m 2 . 

The extended medium simulation was tested against a reference as before. The reference 
mesh was Ax = Ay = As = rjf 16, A z — R/Ns. We generated Ns = 128 independent 
512x512 random phase screens and calculated the field in the receiving plane. The screen 
spacing was then doubled by summing pairs of screens, the field was recalculated and the 
rms intensity error was determined by comparison with the reference field. This procedure 
was repeated to obtain the intensity errors for Az/R — 1 /Ns = 1/64, 1/32, 1/16 and 1/8. 
The rms intensity error is plotted as a function of Az/R in Fig. 4. The estimation error of 
each point is less than the size of the symbol. For a fixed m 2 and pure power-law spectra, 
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Irma has a power-law dependence. It is interesting that I rmt is larger for flatter spectra, 
and also changes more slowly with A z, clearly the high spatial frequency components are 
dominating the intensity errors. For the Kolmogorov exponent a = 5/3, the collection of 
A z errors follow an approximate scaling rule 7/ ms = m 2 Az/47? = ml/ANs-. Although we 
have not advanced a theory for this rule its functional form is not surprising and may be 
useful, particularly when high accuracy is required. When an inner scale is included, the Az 
errors are still lower and, as the screen spacing approaches zero, the error I rm , is linearly 
proportional to Az, in agreement with the results of Spivack and Uscinski 16 . The estimated 
rms intensity m does not change as Az increases, even when I Tms > m. This is very similar 
to the behavior of the transverse sampling error shown in Fig. 1 and likewise implies that 
the error A I due to Az is anticorrelated with the intensity itself. This behavior is consistent 
with the theoretical prediction of Spivack 27 . 

The spatial spectra of intensity and intensity error for two different Az are shown in Fig. 
5. Here we used a pure Kolmogorov spt ctrum and m 2 = 1. It is interesting that the error 
spectra are not white. The error spectra fall at higher wavenumbers such that they never 
exceed the spectrum of intensity by more than a factor of about 2. 

The transverse sampling error for the extended medium simulation was estimated the 
same way as for the thin screen case using 128 screens for each simulation. The results jure 
shown in Fig. 6 for two pure power-law turbulence spectra and a Kolmogorov spectrum 
with an inner scale lo = 5 As. For the pure power-law cases, the best fit power-law for the 
thin-screen result [see Fig. 1] is overplotted showing remarkable agreement. The effect of 
the inner scale is also very similar to the thin screen result. The rms intensity m does not 
change with increasing As until As > 2 s„, as in the thin screen case, indicating that the 
anticorrelation of A 7 and I is similar. 

The spatial spectra of intensity and intensity error due to transverse sampling are shown 
in Fig. 7 for a Kolmogorov spectrum with m 2 = 1. The error spectra are very similar to 
the thin screen spectra shown in Fig. 2, but are not directly comparable because we used a 
smaller m 2 for the extended medium simulation to reduce the computational burden. 

The windowing error is investigated in an analogous manner to the thin screen case. For 
a pure power-law refractive index spectrum and uniform turbulence along the propagation 
path we can use the low wavenumber approximation given by 17 

$/(q) = 87rfc 2 $ n (q) / dzsin 2 [q 2 (7? - z)/(2k)]exp{-D P [j-(R - z)(l +az/R)/(l +a)]}. 

Jo K 

(29) 


The variance error is then 

, a2° +2 r(l +a/2)Dp(r,) 2 r^se 0 RiL f 2 *se 0 R/L f i 

Am = TV1 / dUx / du * / 

flT(l — a/2) Jo Jo Jo 

u~ a ~ 2 sin 2 [tu 2 s 2 /(2r/ 2 )] exp[— u a t a (l -(-a — at)/( 1 + a)]. 


dt x 


(30) 


As with the thin screen case, for fixed Dp(rj) or m 2 the windowing error is a function of 
the parameter 0 o R/L. In the limit of large L the exponential term and the sine term of Eq. 
(30) can be expanded. The leading order result is 


Am 2 a(2 + a)(2ir8rf/ L) 4 °7(a) 
ml 6flT(l — a/2) cos(a7r/4) 


(31) 
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The windowing error, the best fit power-law function, and the theoretical predictions of Eqs. 
(21) and (29) with 8 = 0.65 (the same value as for the thin screen case Fig. 3) are plotted 
in Fig. 8 for the pure Kolmogorov spectrum with m 2 = 1.0. The theoretical scaling law 
agrees with the simulation. However, the statistical fluctuations of the estimates for m 2 are 
much larger than for the thin screen case because fewer realizations could be calculated in 
a reasonable amount of time. 

We can now estimate the total number of mesh points required in a uniform medium 
calculation. If the acceptable error is fixed then the number of points in the three dimensional 
mesh is N p 3 = N p2 R/Az = N p2 Ns so N p 3 a (m 2 ) 4//Q+1 . This very rapid increase of with 
m 2 sets a hard limit on the strength of scattering which can be simulated with this technique. 


7. Point Source in a Uniform Medium 

The case of a point source in an extended medium cannot be analyzed simply by replacing 
the incident plane wave with a point source because of the periodicity implied by the use 
of the DFT. With an incident plane wave the radiation diverges from the axis only due to 
scattering. Thus one period will interfere with another by an overlap of ss 0 O R. However 
in the spherical geometry the radiation is naturally diverging so the overlap from adjacent 
periods is more severe. This problem has been approximated 14,1 using a diverging Gaussian 
beam as the incident field. The average intensity in the observation plane is then a function 
of the radial distance from the beam propagation axis. The statistics of the intensity must 
be calculated in the observation plane using a small region in the center of the simulation 
and normalizing the results by the average intensity. This approach requires careful choice of 
the beam width so the approximation to a spherical wave is good, but the window is not any 
larger than necessary. Another approach, which we have used here, is to cast the split-step 
algorithm in a spherically diverging coordinate system (r,rj x ,7] y ) defined by x = rTj x ,y = rr) y 
for t] Xiy « 1 . 

In this coordinate system the field is written with respect to a monochromatic spherical 


wave 


F{s,r,t) = f,(rf,r)exp[i(kr - ut)]/r. 

(32) 

The parabolic wave equation is then 


df a i (d 2 f s , d 2 f a \ , ., , 

dr = 2*r’( 9, l + d,l) + ' kn(rrh '‘ ^ Z)U 

(33) 


This can be solved in free space (n=0) by Fourier transformation with respect to r/. The 
transform variables are if = s/r and /? = qr. Then 

fl0,r) = /](/?, r 0 )exp[-i/? 2 (l/r 0 - l/r)/(2k)]. (34) 

This expression is used to propagate the field from one screen to the next. The transmittance 
of each screen is defined by its phase spectrum which is given by 

r) = 2 xk 2 / z- 2 <b n (0/z,q z = 0 )dz. (35) 

ir-Ar/2 
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If Ar << r then 

$4,0, r) = 2Kk 2 Ar$ n 0/r,q z = 0 )/r 2 . (36) 

For a pure power-law spectrum ( Iq = 0) and for 0 < a < 2, the strength of scattering 
m 2 , the coherence scale s 0 and the structure function Ds( s) are given by 22 

m 2 b = K 9 (a)D s (rj) (37) 

I< 3 (a) = 2 a r 3 (l + a/2) cos(a7r/4)/r(a + 1) (38) 

D s ( 3) = ( s/s 0 ) a = B 2 (a)k 2 C 2 n Rs a (39) 

B 2 {a) = 2 1 ~ £ T(a)r(l - a/2) sin[(a - l)jr/2]/[(a + l)r(l + a/2)]. (40) 

For the Kolmogorov exponent B 2 ( 5/3) = 1.09289 and K 3 (k>/'3,) — 0.454560. 


The point source simulation was tested against a reference using the same parameters 
and procedure as the plane-wave extended medium case. The intensity errors are plotted 
as a function of Az/R = 1 /Ns in Fig. 9. For a fixed m 2 and pure power-law spectra, the 
error follows a power-law dependence with a similar exponent to that of the plane-wave 
case but the error at a given Az/R is about 4 times higher. For the Kolmogorov exponent 
the collection of intensity errors follows an approximate scaling rule J 2 ma = 4m 2 Az/7? = 
4 ml/ Ns- When an inner scale is included, the errors are smaller and as the screen spacing 
approaches zero, I rm , is linearly proportional to Az , which agrees with the results for the 
plane wave case. 

The spatial spectra of the intensity and the intensity errors for two different values of 
Az are shown in Fig. 10 for the Kolmogorov exponent a = 5/3 and m 2 = 1. The intensity 
error spectra have a higher low-frequency contribution than the plane-wave case. However 
the dominant contribution to the total error variance is still from the higher wavenumber 
region. 

The transverse sampling errors for the point source case are shown in Fig. 11 for two pure 
power-law turbulence spectra and a Kolmogorov spectrum with an inner scale Iq = 5 As. 
The best fit for the thin screen errors shown in Fig. 1 is overplotted. One can see that the fit 
to the pure power-law cases is again remarkably good. Clearly the errors due to transverse 
sampling are well described by the thin screen solution, both for plane and spherical waves. 
The effect of an inner scale is also very similar to that of the thin screen but this case is 
harder to scale. The estimated rms intensity md ec is essentially constant for A s/s 0 < 4. 
The spectral density of the intensity sampling errors are shown in Fig. 12 for a Kolmogorov 
spectrum with m 2 = 1. As with the previous cases, the sampling error spectra are largest 
at high spatial frequencies and display characteristics of aliasing. 

The windowing error is investigated in an analogous manner to the previous two cases. 
For a pure power-law refractive index spectrum and uniform turbulence along the propa- 
gation path, we use the low wavenumber approximation for the intensity spectrum 17,26 and 
calculate the variance error by integrating over the first resolution element. 
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$/(q) = 8irk 2 f dz$ n (qR/z) sin 2 [q 2 2 (i? — z)/{2k)\ exp{— Z?s[q(/2 — z)/k\) (41) 

Jo 


Am 2 = 


2 a+2 a 2 (a + l)r(l + a/2)D s {r } ) 2 [™e 0 R/L ,2*se 0 R/L 


v r(l - a/2) 

ru - ° -2 sin 2 [(l — l/f)u 2 s 2 /(2r/ 2 )] exp[— u"(l — <)"]. 


f2ir6$ 0 R/ L r2x6$ 0 Rf L rl 

I du x I du y dt x 

Jo Jo Jo 


(42) 


For fixed Ds(rj) or m 2 the windowing error is a function of the parameter 0 o R/L. 

In the limit of large L the exponential term and the sin term of Eq. (42) can be expanded. 
The leading order result is 


Am 2 


2r(l + a)(2x(5r / /T) 4 - a /(q) 

ml (a — l)r(l — a/2)T 2 (l + a/2) cos(ax/4) 


(43) 


The windowing error estimates from the simulations, the best fit power-law function, 
and the theoretical predictions of Eqs. (21) and (41) with 8 = 0.75 are plotted in Fig. 8 
for the pure Kolmogorov spectrum with m 2 = 1.0. The theoretical scaling law agrees with 
the simulation at small 0 o R/L. The higher terms of the series expansion for the intensity 
spectrum are required to predict the windowing error for larger 6 0 R/L. 


8. Conclusion 

The FFT split-step algorithm for wave propagation in random media is well-adapted to the 
case of an incident plane wave. It can be applied to the case of an incident spherical wave 
if a diverging spherical coordinate system is used. We find that the error seeding in the two 
cases is essentially the same. The numerical accuracy of the algorithm can be determined 
by the error scaling as a function of transverse sampling (As), longitudinal sampling (As), 
and transverse dimension or window (L). The sampling errors can be studied directly by 
resampling a reference calculation, so one can determine the rms error by comparing the 
resampled calculation with the reference. This allows one to find any of the statistics of 
the error, for example its spectral distribution. The windowing error is more difficult, and 
we have only been able to estimate its effect on the variance by integrating the intensity 
spectrum over a fraction 8 of the spatial frequency bin around zero frequency. Unfortunately 
we cannot create a sample of the windowing error by comparing with a reference field, 
because of the periodicity of the window. 

The transverse sampling error for power-law spectra as a function of 8s/ s 0 is exactly 
the same for the three cases tested, i.e. thin screen and uniform media with plane wave or 
spherical wave incident. It is well described as an aliasing error in the angular spectrum 
so steeper phase spectra have less sampling error than flatter spectra. Phase spectra which 
include an inner scale are particularly steep so they can be simulated with relaxed sampling. 
The spatial spectrum of the sampling error is characterized by maximum contribution at 
the high wavenumbers indicating a small spatial scale and aliasing. 

The longitudinal sampling error is not properly described as an aliasing error although 
it is similar. It is clear the these errors are primarily due to high frequency components 
because steeper power-law spectra have smaller errors as do spectra which include an inner 
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scale. The scaling with A z is similar for plane wave and point source simulations, but the 
point source errors are 4 times higher. The spatial spectrum of the longitudinal sampling 
error is also concentrated at high wave-number but displays a high-frequency cutoff. The 
point-source error spectrum has a larger low-frequency component than the plane- wave case. 

A theoretical expression for the A z scaling of the field error has been given by Spivack 
and Uscinski 16 when the random medium has finite derivatives. The A z error determined 
from the simulations for the plane wave and point source agrees with this scaling when an 
inner scale is included. For a pure power-law spectrum we obtained a weaker scaling law, 
presuemably because the finite derivative constraint was violated. Spivack 27 also calculated 
the ( A 2) 2 dependence for the error of higher moments. This error scaling is difficult to verify 
because the estimation error of the higher moments is large 28,29 . Thus many independent 
samples are required to produce accurate estimates and we were unable to confirm this 
scaling for the intensity variance or any higher moments. 

A similar problem occurred in attempting to verify the theoretical model for the win- 
dowing errors in the intensity variance. For the thin screen simulation, shown in Fig. 3, 
it was easy to calculate many realizations of the field and the windowing error was accu- 
rately determined. For extended medium simulations, which are shown in Fig. 8, the errors 
are considerably larger. However, the windowing error does follow the predictions of our 
theoretical model. 

When an inner scale is included the error scaling is not a simple power-law and the errors 
are less than for the pure power-law case because the high spatial frequency fluctuations of 
the medium and field are reduced. The error analysis presented here can be extended to any 
spectral model for the turbulence and finite beam propagation. It would be very valuable if 
the error in estimating the moments could be expressed in terms of the transverse sampling 
error and the A z error since both of these errors can be accurately estimated with few 
calculations. 
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FIGURES 


Fig. 1. Intensity error I rms due to transverse sampling A s for a plane- wave incident on a 
thin-screen. The phase spectrum for the calculations in the top panel is a pure power-law with 
exponent a = 1 and that in the middle panel has a = 5/3. The phase spectrum used in the lower 
panel is an atmospheric spectrum Eq. (2) with inner scale Iq = 5As. The turbulence level for 
these spectra are described by the Born variance with zero inner scale Eq. (17). These are: m 2 = 
0.1 (o), 1.0 (□), and 10 (O). The best fit power-law is plotted as a straight line on the upper two 
panels. The estimated rms intensity mj ec for = 10 is shown as a dashed line. 

Fig. 2. Spatial spectra of intensity (•) and intensity error due to transverse sampling (open 
symbols) for a plane- wave phase-screen simulation with ( N = 2048, As = tj/ 64). The spectrum 
was a pure Kolmogorov power-law with m£ = 10. The intensity error spectra were calculated for 
As = Tj/ 32 (o) and As = tj/ 16 (□). 

Fig. 3. The intensity variance m 2 (top panel) and the windowing error I rmt (bottom panel) for 
a plane-wave phase-screen simulation as a function of 9 0 R/ L. The phase spectrum was a pure 
Kolmogorov power-law with m 2 = l. The results of the simulation are marked with a (•). The best 
fit power-law is drawn as a solid line and the theoretical model Eqs. (21) and (22) for 6 = 0.65 is 
drawn as a dashed line. 

Fig. 4. The normalized intensity error due to longitudinal sampling A z for a plane- wave extended 
medium simulation. The phase spectra were pure power-law spectra with a = 1 (top) and 5/3 
(middle) and an atmospheric spectrum Eq. (2) with Iq = 5A s (bottom). The turbulence levels 
(for zero inner scale) are marked mj* = 0.1 (o), 1.0 (□), and 10 (O). The normalized rms intensity 
m/mb for m 2 = 10 is shown as a dashed line. 

Fig. 5. Spatial spectra of the intensity (•) and intensity error due to longitudinal sam- 
pling A z (open symbols) for a plane-wave extended medium simulation with ( N = 512, 
As = r//16, Ns = 128). The phase spectrum was a pure Kolmogorov power-law with Born 
variance ra 2 = l. The error spectra were calculated for Ns = 64 screens (o) and 32 screens (□). 

Fig. 6. Intensity error I rm9 due to transverse sampling As for a plane-wave extended medium 
simulation with (N = 512, As = r//16, Ns = 128). The phase spectra were a pure power-law 
with a = 1 (top panel); a = 5/3 (middle) and an atmospheric spectrum Eq. (2) with inner scale 
Iq = 5As (bottom). The turbulence level for these spectra are described by the Born variance with 
zero inner scale Eq. (25). These are: m 2 = 0.1 (o), 1.0 (O), and 10 (O). The best fit power-law 
from the thin screen simulations (Fig. 1) is plotted as a straight line on the upper two panels. The 
estimated rms intensity mj cc for m 2 = 10 is shown as dashed line. 
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Fig. 7. ' Spatial spectra of intensity (•) and intensity error due to transverse sampling (open 
symbols) for a plane-wave extended medium simulation with ( N = 512, As = ry/ 16, Ns = 128). 
The phase spectrum was a pure Kolmogorov power-law with Born variance m£=l. The intensity 
error spectra were calculated for As = rj/8 (o) and As = rj/4 (□). 


Fig. 8. The normalized variance m 2 for plane-wave and point-source extended medium simula- 
tions as a function of 9 0 Rj L. The phase spectrum was a pure Kolmogorov power-law with m 2 = l. 
The best fit power-law models are drawn as solid lines and the theoretical expression are drawn as 
dashed lines. For the plane wave case, the theoretical model is given by Eqs. (21) and (30) with 
6 = 0.65. For the point source case, the theoretical model is given by Eqs. (21) and (42) with 
6 = 0.75. 

Fig. 9. The normalized intensity error due to longitudinal sampling A z for a point-source ex- 
tended medium simulation. The phase spectra were a pure power-law with a = 1 (top) and 5/3 
(middle) and an atmospheric spectrum Eq. (2) with lo = 5As (bottom). The turbulence levels 
(for zero inner scale) are m 2 = 0.1 (o), 1.0 (□), and 10 (O). The estimated rms intensity m/m^ 
for m 2 = 10 is shown as a dashed line. 

Fig. 10. Spatial spectra of the intensity (♦) and intensity error caused by longitudinal sam- 
pling Az (open symbols) for a point-source extended medium simulation with (IV = 512, 
A s = r//16, Ns = 128). The phase spectrum was a pure Kolmogorov power-law with Born 
variance m 2 = l. The error spectra were calculated for Ns = 64 screens (o) and 32 screens (□). 

Fig. 11. Intensity error I rms due to transverse sampling As for a point-source simulation with 
(N = 512, As = rf/l6,Ns = 128). The spectra were a pure power-law spectrum with a = 1 
(top panel); a = 5/3 (middle) and an atmospheric spectrum Eq. (2) with lo = 5As (bottom). 
The turbulence level for these spectra are described by the Born variance (with zero inner scale 
Eq. (37)). These are: mj = 0.1 (o), 1.0 (□), and 10 (O). The best fit power-law from the thin 
screen simulation (Fig. 1) is plotted as a straight line on the upper two panels. The estimated rms 
intensity for m 2 = 10 is shown as a dashed line. 


Fig. 12. Spatial spectra of intensity (•) and intensity error due to transverse sampling (open 
symbols) for a point-source extended medium simulation with (N = 512, As = rj/ 16, Ns = 128). 
The phase spectrum was a pure Kolmogorov power-law with Born variance m 2 = l. The intensity 
error spectra were calculated for As = rj/8 (o) and As = r^/4 (□). 


19 



Figur 


1.0 

0.1 

(0 

E 

k . 

0.01 

0.001 

1.0 

0.1 

<0 

E 

hm 

0.01 

0.001 

1.0 

0.1 

(0 

E 

k. 

0.01 

0.001 




Intensity Spectra 


Figure 2 


Thin Screen Simulation 



Spatial Frequency q 



Am 


Figure 3 




Figure 4 


Plane Wave Extended Medium Simulation 



0.01 0.02 0.04 0.1 0.2 

Az /R = 1/N S 



Figure 5 


Plane Wave Simulation 




rms I rms 

rms 


Figure 6 


Plane Wave Decimation Error 


1.0 

0.1 

0.01 

0.001 

1.0 

0.1 

0.01 

0.001 

1.0 

0.1 

0.01 

0.001 

0 . 




Intensity Spectra 


Figure 7 


Plane Wave Simulation 



Spatial Frequency q 



Figure 8 




rms b rms b 'rms 


Figure 9 


Point Source Simulation 



Az /R = 1/N S 



Intensity Spectra 


Figure 10 


Point Source Simulation 




rms rms rms 


Figure 11 


Point Source Decimation Error 

i 

0 

o.c 

I 

0 

o.c 

I 

0 

o.c 







Intensity Spectra 


Figure 12 


Point Source Simulation 



0.01 0.1 1.0 10.0 

Spatial Frequency q 



